1 Variáveis Pré-Selecionadas

Variável	Descrição
`PropBaixaRenda`	Variável resposta binária (0 ou 1), derivada da “Proporção de pessoas com renda menor de ¼ do salário mínimo”.
`TxAnalfa`	Taxa de analfabetismo do município (%)
`TxDesemp`	Taxa de desemprego do município (%)
`PopRes`	População residente do município (Contagem)
`Saneat`	Quantidade de moradores com abastecimento de água (Contagem)
`ObtOcorr`	Óbitos ocorridos no município

2 Preparação dos Dados

2.1 Leitura

url <- "https://gist.githubusercontent.com/jonhsp/5cf24aa66c3fb173fea117bf4cb1dbc7/raw/29b6a4bcb197ab9be16663f8c73edcfd6b3c51f3/MunicipiosRS.xlsx%2520-%2520Quartis_1_e_3.csv"
df <- read_csv(url)

df %<>%
  select(-'Prop<025SM') %>%
  # Isolando o nome do municipio (remove codigo e espaco inicial)
  mutate(MUNICIPIO = str_remove(MUNICIPIO, "^\\d+\\s")) %>%
  mutate(MUNICIPIO = str_to_lower(MUNICIPIO)) %>%
  # Tratando hifen como zero e convertendo para numerico
  mutate(across(
    .cols = c(TxAnalfa, TxDesemp, PopRes, ObtOcorr, Saneat),
    .fns = ~ {
      .x %>%
        str_replace("^-$", "0") %>%
        str_replace_all("\\.", "") %>%
        str_replace_all(",", ".") %>%
        as.numeric()
    }
  )) %>%
  # Criando a variável binária como fator
  mutate(y_bin = factor(PropBaixaRenda, 
                        levels = c(0, 1), 
                        labels = c("Baixa Proporção", "Alta Proporção"))) %>%
  # Renomeando as colunas para minúsculo
  rename_all(str_to_lower) %>%
  # Selecionando variaveis
  select(municipio, y_bin, txanalfa, txdesemp, popres, obtocorr, saneat) 

summary(df)

##   municipio                     y_bin        txanalfa         txdesemp      
##  Length:249         Baixa Proporção:124   Min.   : 0.900   Min.   : 0.0000  
##  Class :character   Alta Proporção :125   1st Qu.: 3.200   1st Qu.: 0.8525  
##  Mode  :character                         Median : 5.850   Median : 1.7150  
##                                           Mean   : 6.645   Mean   : 2.1223  
##                                           3rd Qu.: 9.700   3rd Qu.: 2.9000  
##                                           Max.   :19.100   Max.   :10.9800  
##                                           NA's   :1        NA's   :1        
##      popres          obtocorr           saneat      
##  Min.   :     0   Min.   :   0.00   Min.   :     0  
##  1st Qu.:  2669   1st Qu.:   6.00   1st Qu.:  2656  
##  Median :  4171   Median :  13.00   Median :  4162  
##  Mean   : 10238   Mean   :  56.28   Mean   : 10183  
##  3rd Qu.:  8771   3rd Qu.:  49.00   3rd Qu.:  8747  
##  Max.   :435564   Max.   :2889.00   Max.   :433266  
##

3 Distribuição das Variáveis

df %>%
  # Seleciona apenas colunas numéricas
  select(where(is.numeric)) %>%
  # Pivota os dados para o formato longo (ideal para facet_wrap)
  pivot_longer(cols = everything(), names_to = "Variavel", values_to = "Valor") %>%
  
  ggplot(aes(x = Valor)) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_density(color = "red", size = 1) +
  # Cria um mini-gráfico para cada variável
  facet_wrap(~ Variavel, scales = "free") +
  labs(
    title = "Distribuição das Variáveis Preditoras",
    x = "Valor da Variável",
    y = "Densidade"
  ) +
  theme_light()

4 Transformações

df_transf <- df %>%
        # Remove municípios com população zero (evita log(0) = -Inf)
        filter(popres > 0) %>%
        # Remoção de NAs
        #drop_na() %>%
        # log populacao
        mutate(log_popres = log(popres)) %>% 
        # Obtos por 100.000 habitantes
        mutate(obt_100k = obtocorr / popres * 100000) %>%
        # Saneamento por 100.000 habitantes
        mutate(saneat_100k = saneat / popres * 100000) %>%
        # Remoção das variáveis originais
        select(-c(popres, obtocorr, saneat))

5 Distribuição das Variáveis Transformadas

5.1 Visualização individual

df_transf %>%
  # Seleciona apenas colunas numéricas
  select(where(is.numeric)) %>%
  # Pivota os dados para o formato longo (ideal para facet_wrap)
  pivot_longer(cols = everything(), names_to = "Variavel", values_to = "Valor") %>%
  
  ggplot(aes(x = Valor)) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_density(color = "red", size = 1) +
  # Cria um mini-gráfico para cada variável
  facet_wrap(~ Variavel, scales = "free") +
  labs(
    title = "Distribuição das Variáveis Transformadas",
    x = "Valor da Variável",
    y = "Densidade"
  ) +
  theme_light()

5.2 Visualização por grupo (Baixa vs Alta Proporção)

# Visualização das variáveis numéricas por nível da variável resposta

df_transf %>%
  # 1. Seleciona numéricas E a variável resposta (y_bin)
  select(y_bin, where(is.numeric)) %>%
  
  # 2. Pivota apenas as numéricas, mantendo o y_bin fixo (cols = -y_bin)
  pivot_longer(cols = -y_bin, names_to = "Variavel", values_to = "Valor") %>%
  
  # 3. Adiciona 'fill = y_bin' para separar as cores por grupo
  ggplot(aes(x = Valor, fill = y_bin, color = y_bin)) +
  
  # position = "identity" permite que os histogramas se sobreponham (transparencia ajuda)
  geom_histogram(aes(y = ..density..), bins = 20, alpha = 0.4, position = "identity") +
  
  # Linha de densidade também separada por cor
  geom_density(alpha = 0.2, size = 1) +
  
  facet_wrap(~ Variavel, scales = "free") +
  
  labs(
    title = "Distribuição das Variáveis por Grupo (Baixa vs Alta Proporção)",
    x = "Valor da Variável",
    y = "Densidade",
    fill = "Grupo",
    color = "Grupo"
  ) +
  theme_light() +
  # Move a legenda para baixo para dar mais espaço aos gráficos
  theme(legend.position = "bottom")

# Visualização das variáveis numéricas por nível da variável resposta (Boxplot)

df_transf %>%
  # 1. Seleciona numéricas E a variável resposta (y_bin)
  select(y_bin, where(is.numeric)) %>%
  
  # 2. Pivota apenas as numéricas, mantendo o y_bin fixo
  pivot_longer(cols = -y_bin, names_to = "Variavel", values_to = "Valor") %>%
  
  # 3. MUDANÇA: Eixo X é o Grupo (y_bin), Eixo Y é o Valor numérico
  ggplot(aes(x = y_bin, y = Valor, fill = y_bin)) +
  
  # Geometria de Boxplot
  geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 1) +
  
  # Escalas livres no eixo Y para acomodar magnitudes diferentes
  facet_wrap(~ Variavel, scales = "free_y") +
  
  labs(
    title = "Boxplots das Variáveis por Grupo",
    x = "Grupo",
    y = "Valor da Variável",
    fill = "Grupo"
  ) +
  theme_light() +
  theme(legend.position = "bottom")

5.3 Relação entre Analfabetismo e Desemprego

# Gráfico de Dispersão: Analfabetismo vs Desemprego por Grupo

df_transf %>%
  ggplot(aes(x = txanalfa, y = txdesemp, color = y_bin)) +
  
  # Adiciona os pontos (dispersão)
  geom_point(alpha = 0.6, size = 2) +
  
  labs(
    title = "Relação entre Analfabetismo e Desemprego",
    subtitle = "Segmentado pela Proporção de Baixa Renda (Y)",
    x = "Taxa de Analfabetismo (%)",
    y = "Taxa de Desemprego (%)",
    color = "Situação (Y)"
  ) +
  
  theme_light() +
  theme(legend.position = "bottom")

6 Ajustes de Modelos

Estamos construindo os modelos progressivamente. O objetivo é encontrar o equilíbrio entre simplicidade e capacidade explicativa.

6.1 Modelo Nulo

O modelo nulo consiste na estimação apenas do \(\beta_{0}\) (intercepto), sem incluir os efeitos das variáveis preditoras. Este modelo será ajustado apenas para comparação da deviance residual.

ajuste_nulo <- glm(y_bin ~ 1, data = df_transf, family = binomial(link = "logit"))

summary(ajuste_nulo)

## 
## Call:
## glm(formula = y_bin ~ 1, family = binomial(link = "logit"), data = df_transf)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.000      0.127       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 343.8  on 247  degrees of freedom
## Residual deviance: 343.8  on 247  degrees of freedom
## AIC: 345.8
## 
## Number of Fisher Scoring iterations: 2

6.2 Modelo 1 - Todas as variáveis

Este é o modelo base (saturado), assumindo que todas as relações são lineares.

# Ajuste do modelo apenas com efeitos principais lineares
ajuste_1 <- glm(y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k + obt_100k, 
                data = df_transf, 
                family = binomial(link = "logit"))

# Sumário para ver coeficientes e significância individual
summary(ajuste_1)

## 
## Call:
## glm(formula = y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k + 
##     obt_100k, family = binomial(link = "logit"), data = df_transf)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.213e+01  1.019e+02   0.413   0.6792    
## txanalfa     1.309e+00  2.008e-01   6.517 7.16e-11 ***
## txdesemp     7.262e-01  2.438e-01   2.979   0.0029 ** 
## log_popres   3.227e-01  4.258e-01   0.758   0.4485    
## saneat_100k -5.364e-04  1.015e-03  -0.528   0.5973    
## obt_100k    -1.951e-03  1.231e-03  -1.585   0.1129    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 343.80  on 247  degrees of freedom
## Residual deviance:  89.62  on 242  degrees of freedom
## AIC: 101.62
## 
## Number of Fisher Scoring iterations: 7

# Análise de Deviance (ANOVA) para ver a contribuição de cada variável ao ser adicionada
# O teste "Chisq" é o adequado para GLMs binomiais
anova(ajuste_1, test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: y_bin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          247     343.80              
## txanalfa     1  230.449       246     113.35 < 2.2e-16 ***
## txdesemp     1   20.541       245      92.81 5.838e-06 ***
## log_popres   1    0.037       244      92.77   0.84706    
## saneat_100k  1    0.140       243      92.63   0.70779    
## obt_100k     1    3.014       242      89.62   0.08257 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.2.1 Diagnóstico do Modelo 1

A análise do modelo inicial (apenas efeitos principais lineares) revela os seguintes pontos:

Qualidade do Ajuste:
- AIC: 101.62 (Base de comparação).
- Deviance Residual: Reduziu de 343.80 (Nulo) para 89.62, indicando que o conjunto de variáveis explica uma grande parcela da variabilidade dos dados.
Significância dos Preditores:
- txanalfa (Analfabetismo): Altamente significativa (\(p < 0.001\)). O coeficiente positivo (\(1.31\)) confirma que quanto maior o analfabetismo, maior a chance de o município pertencer ao grupo de alta pobreza.
- txdesemp (Desemprego): Significativa (\(p = 0.003\)). Coeficiente positivo (\(0.73\)), indicando que o desemprego também impulsiona a pobreza.
- log_popres (População): Não significativa (\(p = 0.45\)). Isso sugere que, em uma relação puramente linear, o tamanho do município não distingue bem os grupos de renda. Isso reforça a necessidade de testar termos quadráticos.
- saneat_100k e obt_100k: Não apresentaram significância estatística neste modelo (\(p > 0.05\)).

Conclusão: O modelo captura bem os efeitos de educação e trabalho, mas falha em capturar a influência da demografia (população) e infraestrutura (saneamento) em sua forma linear.

6.3 Modelo 2 - Termos Quadráticos (Não-Linearidade)

Com base na análise exploratória, adicionamos termos quadráticos para População e Desemprego, mantendo as demais variáveis lineares. O objetivo é testar se a relação com a pobreza é curvilínea.

# Ajuste do modelo com polinômios de 2º grau para as variáveis identificadas
ajuste_2 <- glm(y_bin ~ txanalfa + saneat_100k + obt_100k +
                poly(txdesemp, 2) + 
                poly(log_popres, 2),
                data = df_transf, 
                family = binomial(link = "logit"))

# Sumário para verificar a significância dos termos quadráticos
summary(ajuste_2)

## 
## Call:
## glm(formula = y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 
##     2) + poly(log_popres, 2), family = binomial(link = "logit"), 
##     data = df_transf)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.809e+01  1.023e+02   0.275  0.78369    
## txanalfa              1.332e+00  2.108e-01   6.321  2.6e-10 ***
## saneat_100k          -3.436e-04  1.028e-03  -0.334  0.73817    
## obt_100k             -3.396e-03  1.667e-03  -2.037  0.04170 *  
## poly(txdesemp, 2)1    4.882e+01  1.689e+01   2.891  0.00384 ** 
## poly(txdesemp, 2)2    3.683e+01  1.667e+01   2.209  0.02716 *  
## poly(log_popres, 2)1  1.091e+00  8.316e+00   0.131  0.89559    
## poly(log_popres, 2)2 -2.059e+01  9.309e+00  -2.211  0.02700 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 343.801  on 247  degrees of freedom
## Residual deviance:  80.534  on 240  degrees of freedom
## AIC: 96.534
## 
## Number of Fisher Scoring iterations: 8

# Teste de Deviance (ANOVA) para comparar se o Modelo 2 é superior ao Modelo 1
# H0: Os modelos são iguais (a complexidade extra não ajuda)
# H1: O Modelo 2 é melhor
anova(ajuste_nulo, ajuste_1, ajuste_2, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y_bin ~ 1
## Model 2: y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k + obt_100k
## Model 3: y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) + 
##     poly(log_popres, 2)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)    
## 1       247     343.80                         
## 2       242      89.62  5  254.181  < 2e-16 ***
## 3       240      80.53  2    9.086  0.01064 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3.1 Diagnóstico do Modelo 2 (Termos Quadráticos)

Os resultados obtidos a partir do ajuste do modelo com termos quadráticos indicam uma superioridade clara sobre o modelo linear simples (Modelo 1). A análise detalhada segue abaixo:

Qualidade do Ajuste:
- AIC: 96.534 (Base de comparação).
- Deviance Residual: Reduziu de 89.62 para 80.534, indicando que a inclusão da não linearidade na relação das variáveis preditoras ajuda a explicar melhor os dados, mesmo sendo penalizado pela complexidade adicional dos termos quadráticos..
Teste de Razão de Verossimilhança (ANOVA): A comparação formal entre os modelos resultou em um p-valor de 0.0106 (Chi-quadrado). Como \(p < 0.05\), rejeitamos a hipótese nula de que os modelos são equivalentes. Isso prova estatisticamente que a inclusão da curvatura para Desemprego e População adicionou poder explicativo real ao modelo, e não apenas ruído.
Análise dos Coeficientes:
- Não-Linearidade Confirmada: Os termos quadráticos para Desemprego (poly(txdesemp, 2)2, p=0.027) e População (poly(log_popres, 2)2, p=0.027) foram estatisticamente significativos. Isso sugere que a relação dessas variáveis com a pobreza não é constante, variando conforme a magnitude dos indicadores (formato de curva).
- Preditores Fortes: A Taxa de Analfabetismo (txanalfa) permanece como o preditor mais robusto do modelo (\(p < 0.001\)), mantendo seu efeito positivo forte.
- Variáveis Redundantes: A variável de Saneamento (saneat_100k) apresentou p-valor alto (0.73), sugerindo que, na presença das outras variáveis (especialmente Analfabetismo), ela não adiciona informação nova relevante para a classificação.
- Óbitos: Curiosamente, a taxa de óbitos (obt_100k) tornou-se significativa (\(p=0.04\)) neste modelo ajustado.

Conclusão: O Modelo 2 será adotado como base para a seleção final. O próximo passo será refinar este modelo aplicando o método Stepwise para remover automaticamente os preditores não significativos (como o Saneamento) e obter o modelo final parcimonioso.

7 Seleção do Modelo Final e Diagnóstico

Nesta etapa, refinamos o modelo quadrático removendo variáveis redundantes pelo critério de Akaike (AIC) e realizamos os diagnósticos de qualidade do ajuste.

7.1 Seleção de Variáveis (Stepwise)

Utilizamos o método stepwise (direção “both”) para remover preditores não significativos e chegar ao modelo mais parcimonioso possível.

# Aplica o algoritmo Stepwise no Modelo 2 (Quadrático)
# O R testará remover/adicionar variáveis para baixar o AIC
modelo_step <- step(ajuste_2, direction = "both", trace = 1)

## Start:  AIC=96.53
## y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) + 
##     poly(log_popres, 2)

##                       Df Deviance     AIC
## - saneat_100k          1   80.645  94.645
## <none>                     80.534  96.534
## - poly(log_popres, 2)  2   86.426  98.426
## - obt_100k             1   86.482 100.482
## - poly(txdesemp, 2)    2  101.088 113.088
## - txanalfa             1  276.631 290.631

## 
## Step:  AIC=94.65
## y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) + poly(log_popres, 
##     2)

##                       Df Deviance     AIC
## <none>                     80.645  94.645
## + saneat_100k          1   80.534  96.534
## - poly(log_popres, 2)  2   86.943  96.943
## - obt_100k             1   86.579  98.579
## - poly(txdesemp, 2)    2  101.596 111.596
## - txanalfa             1  277.305 289.305

# Exibe o resumo do modelo vencedor
summary(modelo_step)

## 
## Call:
## glm(formula = y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) + 
##     poly(log_popres, 2), family = binomial(link = "logit"), data = df_transf)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -6.124315   1.168673  -5.240 1.60e-07 ***
## txanalfa               1.334196   0.210857   6.327 2.49e-10 ***
## obt_100k              -0.003382   0.001666  -2.030  0.04237 *  
## poly(txdesemp, 2)1    49.542674  16.900795   2.931  0.00337 ** 
## poly(txdesemp, 2)2    37.308868  16.737594   2.229  0.02581 *  
## poly(log_popres, 2)1   1.500561   8.255232   0.182  0.85576    
## poly(log_popres, 2)2 -20.826304   9.307757  -2.238  0.02525 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 343.801  on 247  degrees of freedom
## Residual deviance:  80.645  on 241  degrees of freedom
## AIC: 94.645
## 
## Number of Fisher Scoring iterations: 8

O modelo capta a não linearidade dignificativa de log populção e da taxa de desemprego, porém também acrescenta a taxa de obtos ao nível de 5% de significancia, porém a mesma não seria incluida ao nível de 1% de significancia. O \(\beta_{obtos}\) foi estimado em -0.003, o que indica que a taxa de obtos tem um efeito negativo na probabilidade de alta taxa proporção de pobreza. Esse resultado calsa estranheza, pois indica que municípios mais violetos tendem a ter menores taxas de pobreza, contudo, a significância desta variável deve estar relacionada com o efeito de outros fatores, externos ao modelo, que parcialmente são representandos pela quantidade de obtos.

Visando a simplicidade do modelo e a interpretabilidade dos resultados, a taxa de obtos será removida, sem grandes perdas ao modelo, visto que (AIC) < 4

modelo_final <- update(modelo_step, . ~ . - obt_100k)

# Exibe o resumo do modelo vencedor
cat("----Summary do Modelo Final----")

## ----Summary do Modelo Final----

summary(modelo_final)

## 
## Call:
## glm(formula = y_bin ~ txanalfa + poly(txdesemp, 2) + poly(log_popres, 
##     2), family = binomial(link = "logit"), data = df_transf)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -7.1293     1.1086  -6.431 1.27e-10 ***
## txanalfa               1.2417     0.1895   6.553 5.66e-11 ***
## poly(txdesemp, 2)1    42.5755    15.1374   2.813  0.00491 ** 
## poly(txdesemp, 2)2    29.5031    14.6302   2.017  0.04374 *  
## poly(log_popres, 2)1  -5.4882     8.1155  -0.676  0.49888    
## poly(log_popres, 2)2 -13.3643     8.9366  -1.495  0.13479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 343.801  on 247  degrees of freedom
## Residual deviance:  86.579  on 242  degrees of freedom
## AIC: 98.579
## 
## Number of Fisher Scoring iterations: 8

cat("----Anova entre os Modelos----")

## ----Anova entre os Modelos----

anova(modelo_final,modelo_step, ajuste_2, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y_bin ~ txanalfa + poly(txdesemp, 2) + poly(log_popres, 2)
## Model 2: y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) + poly(log_popres, 
##     2)
## Model 3: y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) + 
##     poly(log_popres, 2)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       242     86.579                       
## 2       241     80.645  1   5.9332  0.01486 *
## 3       240     80.534  1   0.1115  0.73850  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("---- AIC ----")

## ---- AIC ----

AIC(ajuste_nulo, ajuste_1, ajuste_2, modelo_step, modelo_final)

##              df       AIC
## ajuste_nulo   1 345.80100
## ajuste_1      6 101.62016
## ajuste_2      8  96.53401
## modelo_step   7  94.64546
## modelo_final  6  98.57871

7.2 Interpretação dos Efeitos (Odds Ratios)

A tabela abaixo apresenta a Razão de Chances (Odds Ratios) para as variáveis do modelo final.

# Extração dos coeficientes e cálculo do OR (exp(beta))
tabela_or <- exp(cbind(OR = coef(modelo_final), confint(modelo_final)))

# Exibição formatada
knitr::kable(tabela_or, digits = 3, caption = "Razão de Chances (Odds Ratios) e IC 95%")

Razão de Chances (Odds Ratios) e IC 95%
	OR	2.5 %	97.5 %
(Intercept)	1.000000e-03	0.000	5.000000e-03
txanalfa	3.462000e+00	2.508	5.334000e+00
poly(txdesemp, 2)1	3.092593e+18	55970162.875	1.107302e+33
poly(txdesemp, 2)2	6.501656e+12	143.121	4.430122e+26
poly(log_popres, 2)1	4.000000e-03	0.000	2.614928e+04
poly(log_popres, 2)2	0.000000e+00	0.000	1.210000e+01

8 Visualizações

8.1 Visualização ( Variáveis fixadas na média)

# 1. Efeito do Analfabetismo (Linear)
# Mostra a probabilidade predita de Alta Pobreza conforme a Taxa de Analfabetismo aumenta
plot_analfa <- ggpredict(modelo_final, terms = "txanalfa") %>%
  plot() +
  labs(
    title = "Efeito do Analfabetismo",
    x = "Taxa de Analfabetismo (%)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# 2. Efeito do Desemprego (Quadrático)
# Mostra a curva não-linear capturada pelo polinômio
plot_desemp <- ggpredict(modelo_final, terms = "txdesemp [all]") %>%
  plot() +
  labs(
    title = "Efeito do Desemprego (Não-linear)",
    x = "Taxa de Desemprego (%)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# 3. Efeito da População (Quadrático)
# Plota a curvatura ajustada para o Log da População
plot_pop <- ggpredict(modelo_final, terms = "log_popres [all]") %>%
  plot() +
  labs(
    title = "Efeito do Tamanho da População",
    x = "Log(População)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# Exibir os gráficos juntos usando gridExtra
# Organiza em 2 colunas: Analfabetismo e Desemprego na linha 1, População centralizada na linha 2
grid.arrange(plot_analfa, plot_desemp, plot_pop, ncol = 3, 
             layout_matrix = rbind(c(1, 2), c(3, 3)))

8.2 Visualização ( Variáveis fixadas na moda)

# 1. Efeito do Analfabetismo (Linear)
# Mostra a probabilidade predita de Alta Pobreza conforme a Taxa de Analfabetismo aumenta
plot_analfa <- ggpredict(modelo_final, terms = "txanalfa", typical = 'mode') %>%
  plot() +
  labs(
    title = "Efeito do Analfabetismo",
    x = "Taxa de Analfabetismo (%)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# 2. Efeito do Desemprego (Quadrático)
# Mostra a curva não-linear capturada pelo polinômio
plot_desemp <- ggpredict(modelo_final, terms = "txdesemp [all]", typical = 'mode') %>%
  plot() +
  labs(
    title = "Efeito do Desemprego (Não-linear)",
    x = "Taxa de Desemprego (%)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# 3. Efeito da População (Quadrático)
# Plota a curvatura ajustada para o Log da População
plot_pop <- ggpredict(modelo_final, terms = "log_popres [all]", typical = 'mode') %>%
  plot() +
  labs(
    title = "Efeito do Tamanho da População",
    x = "Log(População)", 
    y = "Probabilidade de Alta Pobreza"
  ) +
  theme_light()

# Exibir os gráficos juntos usando gridExtra
# Organiza em 2 colunas: Analfabetismo e Desemprego na linha 1, População centralizada na linha 2
grid.arrange(plot_analfa, plot_desemp, plot_pop, ncol = 3, 
             layout_matrix = rbind(c(1, 2), c(3, 3)))

8.3 Visualização 3D das Probabilidades

# 1. Preparar os dados com as predições
# Usamos os dados reais para ver onde os municípios se posicionam no espaço 3D
df_3d <- df_transf %>%
  mutate(
    prob_predita = predict(modelo_final, type = "response"),
    # Criamos um texto formatado para quando passar o mouse
    hover_text = paste(
      "<b>Município:</b>", municipio,
      "<br><b>Analfabetismo:</b>", txanalfa, "%",
      "<br><b>Desemprego:</b>", txdesemp, "%",
      "<br><b>População:</b>", exp(log_popres) %>% round(0), # Revertendo o log para ler fácil
      "<br><b>Prob. Alta Pobreza:</b>", scales::percent(prob_predita, accuracy = 0.1)
    )
  )

# 2. Gerar o Gráfico de Dispersão 3D
plot_ly(data = df_3d, 
        x = ~txanalfa, 
        y = ~txdesemp, 
        z = ~log_popres,
        color = ~prob_predita, 
        colors = c("#440154", "#21908C", "#FDE725"), # Escala Viridis (escura -> clara)
        type = "scatter3d", 
        mode = "markers",
        marker = list(
          size = 6,              # Tamanho das bolinhas
          opacity = 0.8,         # Transparência para ver pontos atrás
          line = list(width = 0) # Remove borda para ficar mais limpo
        ),
        text = ~hover_text,
        hoverinfo = "text"
) %>%
  layout(
    title = "Espaço Tridimensional dos Preditores<br><sub>Cor = Probabilidade de Alta Pobreza</sub>",
    scene = list(
      xaxis = list(title = "Analfabetismo (%)"),
      yaxis = list(title = "Desemprego (%)"),
      zaxis = list(title = "Log(População)"),
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.2)) # Posição inicial da câmera
    )
  )

8.4 Visualização 3D de Acertos e Erros

# 1. Preparar os dados de classificação
df_classificacao <- df_transf %>%
  mutate(
    # Probabilidade predita pelo modelo
    prob = predict(modelo_final, type = "response"),
    
    # Classificação baseada no limiar de 50% (0.5)
    # Se prob > 0.5 -> Classifica como "Alta Proporção" (nível 2 do fator)
    # Se prob <= 0.5 -> Classifica como "Baixa Proporção" (nível 1 do fator)
    classe_predita = ifelse(prob > 0.5, "Alta Proporção", "Baixa Proporção"),
    
    # Verificar se acertou ou errou
    # Como y_bin é fator, convertemos para caracter para comparar
    resultado = ifelse(as.character(y_bin) == classe_predita, "Acerto", "Erro"),
    
    # Detalhe do erro (Falso Positivo vs Falso Negativo) para o tooltip
    tipo_erro = case_when(
      resultado == "Acerto" ~ "Classificação Correta",
      classe_predita == "Alta Proporção" ~ "Falso Positivo (Predisse Alta, mas é Baixa)",
      TRUE ~ "Falso Negativo (Predisse Baixa, mas é Alta)"
    ),
    
    # Texto para o mouse
    hover_text = paste(
      "<b>Município:</b>", municipio,
      "<br><b>Analfabetismo:</b>", txanalfa, "%",
      "<br><b>Desemprego:</b>", txdesemp, "%",
      "<br><b>População:</b>", exp(log_popres) %>% round(0), # Revertendo o log para ler fácil
      "<br><b>Prob. Alta Pobreza:</b>", scales::percent(prob, accuracy = 0.1)
    )
  )

# 2. Gerar o Gráfico 3D
plot_ly(data = df_classificacao, 
        x = ~txanalfa, 
        y = ~txdesemp, 
        z = ~log_popres,
        color = ~resultado, 
        colors = c("#1f77b4", "#d62728"), # Azul para Acerto, Vermelho para Erro
        symbol = ~resultado,
        symbols = c("circle", "circle"), # Bolinha para Acerto, X para Erro
        type = "scatter3d", 
        mode = "markers",
        marker = list(
          size = 6,
          opacity = 0.8,
          line = list(width = 1, color = "black") # Borda para destacar
        ),
        text = ~hover_text,
        hoverinfo = "text"
) %>%
  layout(
    title = "Diagnóstico de Classificação 3D<br><sub>Pontos vermelhos indicam onde o modelo errou</sub>",
    scene = list(
      xaxis = list(title = "Analfabetismo (%)"),
      yaxis = list(title = "Desemprego (%)"),
      zaxis = list(title = "Log(População)")
    ),
    legend = list(title = list(text = "Resultado"))
  )

8.5 Visualização da probabilidade conjunta

# Definindo as medianas para fixar as variáveis de controle
med_analfa <- median(df_transf$txanalfa)
med_desemp <- median(df_transf$txdesemp)
med_pop    <- median(df_transf$log_popres)

# ----------------------------------------------------------------------------
# GRÁFICO 1: Analfabetismo vs Desemprego (Fixando População)
# ----------------------------------------------------------------------------
grid1 <- expand.grid(
  txanalfa = seq(min(df_transf$txanalfa), max(df_transf$txanalfa), length.out = 50),
  txdesemp = seq(min(df_transf$txdesemp), max(df_transf$txdesemp), length.out = 50),
  log_popres = med_pop
)
grid1$prob <- predict(modelo_final, newdata = grid1, type = "response")

p1 <- ggplot(grid1, aes(x = txanalfa, y = txdesemp)) +
  geom_tile(aes(fill = prob)) +
  geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
  geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) + # Pontos reais
  scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1)) +
  scale_color_manual(values = c("black", "darkred")) +
  labs(title = "Analfabetismo vs Desemprego", subtitle = "População fixa na mediana") +
  theme_minimal() + theme(legend.position = "none")

# ----------------------------------------------------------------------------
# GRÁFICO 2: Analfabetismo vs Log(População) (Fixando Desemprego)
# ----------------------------------------------------------------------------
grid2 <- expand.grid(
  txanalfa = seq(min(df_transf$txanalfa), max(df_transf$txanalfa), length.out = 50),
  log_popres = seq(min(df_transf$log_popres), max(df_transf$log_popres), length.out = 50),
  txdesemp = med_desemp
)
grid2$prob <- predict(modelo_final, newdata = grid2, type = "response")

p2 <- ggplot(grid2, aes(x = txanalfa, y = log_popres)) +
  geom_tile(aes(fill = prob)) +
  geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
  geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) +
  scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1)) +
  scale_color_manual(values = c("black", "darkred")) +
  labs(title = "Analfabetismo vs População", subtitle = "Desemprego fixo na mediana") +
  theme_minimal() + theme(legend.position = "none")

# ----------------------------------------------------------------------------
# GRÁFICO 3: Desemprego vs Log(População) (Fixando Analfabetismo)
# ----------------------------------------------------------------------------
grid3 <- expand.grid(
  txdesemp = seq(min(df_transf$txdesemp), max(df_transf$txdesemp), length.out = 50),
  log_popres = seq(min(df_transf$log_popres), max(df_transf$log_popres), length.out = 50),
  txanalfa = med_analfa
)
grid3$prob <- predict(modelo_final, newdata = grid3, type = "response")

p3 <- ggplot(grid3, aes(x = txdesemp, y = log_popres)) +
  geom_tile(aes(fill = prob)) +
  geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
  geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) +
  scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1), name="Prob.\nPredita") +
  scale_color_manual(values = c("black", "darkred"), name="Real") +
  labs(title = "Desemprego vs População", subtitle = "Analfabetismo fixo na mediana") +
  theme_minimal() + theme(legend.position = "right")

# ----------------------------------------------------------------------------
# Exibição Conjunta
# ----------------------------------------------------------------------------
# Organiza os 3 plots. O p3 fornece a legenda para todos.
grid.arrange(p1, p2, p3, ncol = 3)

9 Diagnóstico do Ajuste

9.1 Multicolinearidade (VIF)

car::vif(modelo_final)

##                         GVIF Df GVIF^(1/(2*Df))
## txanalfa            1.217938  1        1.103602
## poly(txdesemp, 2)   1.697873  2        1.141501
## poly(log_popres, 2) 1.591663  2        1.123215

10 Teste de Shapiro-Wilk

shapiro.test(residuals(modelo_final))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_final)
## W = 0.84722, p-value = 6.342e-15

O GVIF, bem como o teste de Shapiro Wilk, confirmam que não há evidência de multicolinearidade entre as variáveis preditoras.

11 Resíduos

hnp(modelo_final, resid.type = "deviance", how.many.out = TRUE, paint.out = TRUE,
    main = "Envelope Simulado", xlab = "Percentil Teórico", ylab = "Resíduos Deviance")

## Binomial model

## Total points: 248 
## Points out of envelope: 0 ( 0 %)

plot(residuals(modelo_final))

Trabalho 1: Regressão Logística - Proporção de Baixa Renda (RS)

CE314 - Modelos Lineares Generalizados (Grupo 4)

jonathas Pereira,……